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TECHNICAL MEMORANDUM X-53459 


VARIATIONAL PROBLEMS AND THEIR SOLUTION BY THE ADJOINT METHOD 

SUMMARY 


A variational problem typical of those encountered in flight mechanics 
is posed. The adjoint technique is then developed in a manner to indicate 
its application to general two-point boundary value problems. It is then 
specifically applied to the variational problem. Finally, the practicality 
of the adjoint method is illustrated by solving three typical problems of 
exo-atmospheric flight and one problem involving an ascent through the 
atmosphere to low earth orbit. 


I. INTRODUCTION 


Many flight mechanical problems require that some "cost" criterion 
be extremized. A typical criterion is that the scientific payload of a 
rocket booster be maximized. Besides maximizing the payload, certain 
functional relationships among the "state" variables of the vehicle must 
invariably be satisfied at the end point. The calculus of variations, a 
branch of classical mathematics under development for over 250 years, 
forms the basis for handling optimization problems of this type. Unfor- 
tunately, this mathematical treatment leads to two-point boundary value 
problems which in themselves are important stumbling blocks to a straight- 
forward application of the theory. The adjoint method has been developed 
as a tool to solve the two-point boundary value problem which in turn 
allows the ideas and results of the classical calculus of variations to 
be applied in a great variety of interesting and practical problems. The 
remainder of this text is devoted to formulating the variational problem, 
developing the adjoint method of solution to the split boundary value 
problem, and applying the resultant theory numerically to some typical 
trajectory problems. No claim for mathematical rigor is made. Necessary 
continuity requirements and the existence of derivatives of the required 
order are assumed. The treatments and developments of the adjoint method 
are for the most part formal in nature but their utility and worth have 
been borne out by the solution of a number of practical problems. 



II. THE VARIATIONAL PROBLEM 


Consider a set of ordinary first order differential equations: 


*i = ••• » x m > u i> ••• > u n ) i = 1, , m (2.1) 


= f i (x,u) 


where the defining the "state 11 of the system are the dependent variables 
and the u l5 ... , u n are the control or forcing variables which are implicit 
functions of t, the independent variable, which here is taken as time. A 
solution of this system in an interval t G ^ t ^ tf is given by m functions, 
X£(t), and n functions, uj(t), such that their substitution reduces equa- 
tions (2.1) to identities . The system is said to have n degrees of freedom. 
The state values at t Q , together with t 0 , are termed the initial boundary 
and the state values at tf together with tf are called the terminal boun- 
dary. The variational problem consists of extremizing a functional 


J = 0(xi, 


*m) t f + 



f 0 < x i> 


» ^Sn’ u i’ 


= 0(x) 


C f 



f 0 (x,u) dt 


u n ) dt 


( 2 . 2 ) 


subject to the differential constraints (2.1) and the terminal constraints 

£ ( X l> •••> Xjjj) ^ — 0 Z ~ 1» ...,q (2.3) 

q s m + 1 


where the are given constants. Employing the method of Lagrange, an 
extremal requires that the first variation of the following functional 
vanish: 
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J = 0( Xl , .... + J f o (x x , .... x^.Ui, ..., i^) 

fc o 


+ q i ( f i (x> • V u i’ •••» "n> " *i)J dt + v £ ^(x 1( 

.... x m ) tf J 


(2.4) 

where it is understood that a repeated subscript implies summation over 
its range, the q^(t) are Lagrange multiplier functions, and v. are Lagrange 
multiplier constants. 

The vanishing of 6J leads to the following necessary conditions: 

. _ dH \ . . 

1,1 ' ' aVu,, 

(2.5) 

sir) ■ 0 J * 1 n 

(2.6) 



where 


H(x,u,q) - f Q (x,u) + q i f i (x,u) 

(2.7) 

is the Hamiltonian, 


q ‘ (t£> * © tf ' %) tf 

(2.8) 

£ o (x,u) + - [v, ^ . 

f C f 

(2.9) 

i = 1, . . . , m 
l = 1, . . . , q 
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Notice that H is a constant of the system since 


dH dH A • 

dt bx.J X i Sq. 

i/u,q n i/x,u 


. , SH \ 

<1 i + 5rj u j 

j/x,q J 


Substituting (2*1), (2.5), and (2.6) further yields 


= -q.x. + f .q. + 0 = -q.x. + q.x. = 0. 
dt x 1 i n i n i 1 1 


The necessary conditions given have implicitly assumed that the 
initial state boundary is known. Equations (2.5) and (2.6) are the 
usual Euler-Lagrange equations. Equations (2.8) and (2.9) express the 
transversal ity conditions at the terminal boundary. A solution to the 
variational problem consists of finding functions x(t), u(t) and q(t) 
in the interval t D ^ t ^ tf such that equations (2.1) and (2.5) are 
reduced to identities, equation (2.6) is satisfied and the boundary con- 
ditions, equations (2.3), (2.8) and (2.9) are satisfied. The u(t) are 
implicitly defined by (2.6) which means that solutions of the differential 
system (2.1) and (2.6) are determined by 2m + 2 boundary values of which 
m + 1 values are provided by the assumption that the initial state bound- 
ary is known and the remaining m + 1 values are obtained from the 
m + q + 1 terminal conditions, (2.3), (2.8), and (2.9) by eliminating 
v i> •••> Vq .* Thus, an extremal requires the solution of a two-point 
boundary value problem. The nature of such a problem, as well as the 
fact that (2.1) and (2.5) are normally highly nonlinear, means that an 
analytic solution is rarely possible and a numerical solution must be 
sought. Numerical solutions require that one complete boundary be known. 
The classical technique of solution is to guess the unknown initial bound- 
ary values and to adjust them until the terminal boundary conditions are 
satisfied. This transforms the split boundary value problem into an init- 
ial value problem. The manner of solution described is called the indirect 
method since the control variables are obtained implicitly as functions of 
time and are not specified explicitly and then modified in a manner such 
that (2.3) is satisfied while simultaneously extremizing (2.2). Steepest 
descent techniques are of this latter type which are called direct methods. 


If (2.3) are m in number, then m of (2.8) and (2.9) are superfluous. 
If (2.3) are m + 1 in number, both (2.8) and (2.9) are unnecessary. 
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III. THE ADJOINT SOLUTION TO TWO- POINT BOUNDARY VALUE PROBLEMS 


Consider a differential system of the following type: 


* V' 1 - 


* 2 »> 


a = 1, 


2m 


(3.1) 


V z) - 


The F(x are supposed to be functions of the differential dependent variables 
only and the conments concerning (2.1), except for the control variables, 
are applicable. Let 


5z a (t) = z a (t) ‘ z a (t) ’ 


(3.2) 


where z a (t) is a nominal solution of (3.1), z a (t) is a varied solution 
and the concept of a variation symbolized by 5 is discussed in Appendix I. 
The variations in z(t) obviously induce (to first order) a variation in 
the F a (t) as follows: 


5F a = F a^ z + 5z ^ " F a^ * cteT 5z p* ^ * 1 » •••» 2,11 ( 3 * 3 ) 


Combining (3.1) and (3.3) yields 


j dF 

5 ^a> - £ < 8z *> - 8F a = ^ 5z 3 

P 


or 


57 < 6 V ' § 

P 

where the interchange of d/dt and 6 are justified in Appendix I. 


(3.4) 

The 


5 



equations adjoint to (3.4) are defined as 


df 

p a = • v a,p = 


2m 


(3.5) 


where the p a are the adjoint variables. The usefulness of the adjoint 
variables (also called the influence functions) is indicated by the 
following: 


K bz a> = P, a Sz a + p a (8z a> 

“ ' 3^ Pp 8z a + Pa 3^ 6z p 


= 0 


(3.6) 


where (3.4) has been substituted and the fact that oc and p are dummy sum- 
mation symbols has been used. Since p a bz a is a constant of the system, 
obviously " 

p a * p a 5 z a)- ■ < 3 - 7 > 

o t t 


where t£ is the nominal terminal time. Herein lies the usefulness of the 
adjoint variables; they relate variations at the nominal terminal boundary 
to variations at the initial boundary. Again let boundary conditions of 
the following form be imposed at the initial and terminal boundaries: 


(i) (z) t = o 

i~ 1, . • • y P 

(3.8a) 

O 



(j) (z) t = 0 . 

j = p + 1, ..., 2m - p + 1 

(3.8b) 


6 



Boundary values on the adjoint variables have not been imposed. A 
meaningful interpretation can be attached to the following choice for 
the set of boundary values at the nominal terminal time: 



j = p + 1, .... 2m - p + 1 (3.9a) 

CK = 1 , . . . , 2m 


As shown in Appendix I, at the nominal terminal time, 



+ z 


a 



(3.10) 


where dz a is the total differential change in Zq, at tf. Thus, 





( dz a ' 2 a dt) 





( d2 a - dt) 



(3.11) 


Substituting the relations 


dfl 


(j) _ dft 


(j) 


dz 


a 


0) _ M 


a 

(j) 


a 


5z a a 


into (3.11) results in 


(3.12) 




U <3) - fi <3) 



(3.13) 
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Substituting (3.13) into (3.7) yields 






dt 


-it 


f 


or 


,U> 

a 


Sz 


a 


+ dt 

= dft^ 

_t o 

_t f J 


(3.14a) 


p of the Sz)t Q are eliminated in (3.14a) by using (3.8a). Assuming the 

are specified, (3.14a) represents 2m - p + 1 equations for one 
fc f 


unknown dtf and 2m - p unknowns, Bz.)t • Equations (3.14a) are generated 
x 1 o 

by 2m - p + 1 integrations of the adjoint equations. This number may be 
reduced by one if the nominal terminal time is determined by one of the 

terminal boundary conditions being satisfied, say P + ^(z)£ = 0. 

Then, 


(2m-p+l), . ^(2m p+i) 

d^ (z); - — dz a 

fc f oz a J t, 


=0. i = 1 , ..., 2m 


If the adjoint equations are integrated with boundary conditions, 


a A 


Sn 


(j) 


sr 


(j) /^(2m-p+l) x -i 


Bz o: a (2m " p+1) 


"3£ 


a /Jt 


(3.9b) 


a straightforward calculation shows that (3.14a) reduce to 


Sz 1 = 


j = 1 , . . . , 2m - p 


(3.14b) 



Of course, these last operations are not possible if the designated 
terminal boundary condition cannot be satisfied from the initial boundary. 

The implementation of the adjoint method is an iterative process. 

This means that each successive step is dependent on the previous step. 

The process is said to have converged when the differential equations and 
the terminal boundary conditions are reduced to identities by the solution 
functions. When t G is fixed, a nominal solution z(t) is generated by 
assuming a nominal t f , guessing the missing initial boundary values, and 
integrating (3.1) forward in time until t = t f . Normally the 



and changes in tf and the free initial boundary values are necessary. 
These changes are computed by using the initial conditions of (3.9a) to 
integrate the adjoint equations (3.5) backward in time until t = t D . 
This yields the 


P 


(j) 


( fc o) 


needed for (3.14a). Since (3.14a) is generally the result of a linear 
treatment of a nonlinear system, it is unreasonable to expect it to pre- 
dict the whole correction necessary during the early steps of the con- 
vergence process. Consequently, only fractional portions of the terminal 
boundary condition violations are used in (3.14a) initially; i.e., the 


are replaced by 




2m - p + 1 


where 0 < Cfcs l, As convergence is attained, the are increased 
since the linear approximation is becoming better and better. Using these 
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(3.14a) is used to calculate the new initial boundary as 


new 


z. J = z. I + Sz 

xJ t V 

L o L o 


old 


and the new terminal time as 


jnew . t old + d j £ 



The state equations are again integrated forward in time until t = tf 
Now the 


dft 


u>y 


new 


rnew 

c f 


will have either increased or decreased. 

(1) If they have increased, the changes calculated are too 
large and the linearity assumptions have obviously been 
violated. This is corrected by halving the computed 
changes and reintegrating the state equations forward. 

This secondary iteration is done as many times as required 
to decrease the terminal boundary condition violations. 
When they have finally been reduced, the adjoint equations 
are again integrated backwards according to the previous 
step. 

(2) If they have decreased, the adjoint equations are immedi- 
ately integrated backwards and the previous step repeated. 

This process continues until the terminal boundary conditions are satis- 
fied to within some tolerance. The secondary iteration described above 
is essential since it can reduce the number of times the adjoint equa- 
tions are integrated. This results in a large time saving since these 
equations are usually rather lengthy and involved. 
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IV. THE ADJOINT SOLUTION OF THE VARIATIONAL PROBLEM 


Seeking to limit the scope of the discussion while at the same time 
illustrating the application of the adjoint technique to a problem of 
some interest, a particular type of variational problem is considered 
here: the time minimal problem. The functional to be minimized is, 

thus. 


J = 




(4.1) 


i.e., f 0 (xi» ..., xjujUj, ..., i^) = 1. Further assumptions are that the 
initial boundary and the terminal boundary (except for tf) are known. 

The problem, then, is to move from an initial known state to a terminal 
known state in the least time. In order that the scheme outlined in Sec- 
tion III be applicable to this variational problem, it is necessary that 
the dependence of the differential equations on the control variables be 
removed. This can be done in principle by using (2.6) where the relations 
are made explicit as 


Uj = Uj(x,q). j = 1, ..., n 

The differential equations of this variational problem become 



subject to the boundary conditions 



(4.2) 


(4.3) 


(4.4 


(4.5) 


(4.6) 
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(4.7) 


q 


+ v. 



= o 


1 


+ q fj*. 


u ( x »q) 



f 


= o. 


(4.8) 


Equations (4.3) and (4.4) can be written in the form (3.1) by making the 
identifications 


z a s x i’ 

z a 

•rl 

4-1 

III 

a = 

1 , • . . , m 

z a s q i* 


III 

0Q 

» 

a = 

m + 1 , ...» 2m 


Then, 


a 


= F a (x,q). 


a * 1, . . . , 2m. 


(4.9) 


Boundary conditions (4.6), (4.7), and (4.8) are m + 1 in number if the 
Vi, ..., Vq are eliminated and correspond to (3.8) which are m + 1 in 
number since the initial state boundary has been assumed known. Con- 
sequently, the formulation of Section III is applicable. The adjoint 
equation (3.5) may be conveniently rewritten with indices varying from 
1 to m as 


.x 


Pi 



(4.10a) 



(4.10b) 


where the superscripts have been used to indicate the correspondence 
between the subscripts of (3.5) and the identifications in (4.9). 
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To calculate the partial derivatives in (4.10), the following differentials 
of (4.2), (4.3) and (4.4) are required: 


du. 

J 




(4.11) 


^ f i\ 

df i ' Kl ** 



i,k = 1, 

j - 1, 


m 

n 


(4.12) 


d *i 



+ 


^Vx.u 


d v 


(4.13) 


Substituting (4.11) into (4.12) and (4.13) and rearranging yields 



(4.15) 


By holding x or q constant as required, (4.14) and (4.15) can be used to 
calculate the partial derivatives in (4.10). The results are 


.x 


P 


i 




(4.16) 



(4.17) 
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The only unknowns in the last two equations are 


se) 

i/q 


and 


du. 

dq.y 

l'x 


These partial derivatives may be calculated from (4.5) as follows: 


d f A 

i r (x,u,q) = q ± - 0 

r 7 x 


ah A * r \ 8*> t \ 

dh r ■ sr) d *k + sr) du j + sr) d< > 

K'u,q y x,q J ^k'x.u 


0. (4.18) 


Holding q constant in (4.18) yields 


dh \ dh \ du.\ • t i 

sr) + sr) s£) ■ °> ' " 

i/u,q j 7 x,q t'q j,r = 1, n 


(4.19) 


which represents n equations in the unknowns 


du.\ 

^) q ' 

Consequently, if the array 


3h r\ 

^x„ 


where r is the row index and j is the column index - is nonsingular. 



Holding x constant in (4.18) yields 


With the aid of 
(4.17) become 




(4.21) 


(4.20) and (4.21), the adjoint equations (4.16) and 



(4.22) 



(4.23) 
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Equations (4.3) - (4.8), (4.22) and (4.23) can be expressed more succrntly 
in terms of the Hamiltonian which, for this time minimal problem, is 


H(x,u,q) = 1 + q i f i . 


(4.24) 


These equations can then be written 


. _ c>H 

X i " Sq . . 

1 ' x,u 


(4.25) 


• _ SH 

q i dx. 

1 q , u 


(4.26) 


^H . _ _ 

S7T j - 0 

jx.q 


(4.27) 


- Ej = 0 

Z f 


(4.28) 


q i + v t —) t = 0 

1/ t c 


(4.29) 


H )t f = 0 


(4.30) 


• x _ 

3 . ~ 

’ 1 


<3 2 h ^ 

_ S 2 h ^ 

u,x,q 

B 2 H 

S 2 h \ 

^ x i 

K,u,q V 

5u . ()u 
- J rJ 

x, q , x ^*i ^Vx,q,u 


d 2 H 


S 2 H 


Vu,q,U Vu,q,x 



S 2 H ] 

u,q,X 

ciu . ckl 

L J r J 


S 2 h 


'x, q , x ^ x i ^Vx,q,u 



.q _ 

p l~ 


S g H 

^7^ 


\ 

d a H 

) 

u,x,q 

L J rJ 


-,-1 


S 2 h 


du dq. 

r 1 ' u,q,x 


d 2 H 


d g H 


dq.dx./ du. bx.J 

i Vq,u,x 2 i'u,q. 


d 2 H 

d 2 H \ 

du. du 
J r-J 

x,q,x ^k'u,q,x 


P k* 


(4.32) 


where the nonstandard triple subscript notation indicates, for example, 
that 


o g H 


is to be evaluated by finding the partial derivative of H with respect 
to holding x and u constant and then finding the partial derivative 
of this result with respect to X£ by holding u and q constant. Equa- 
tions (4.24) - (4.32) are convenient when an analytic derivative cal- 
culator such as IBM's FORMA C system is available, since all differential 
equations are essentially generated by the Hamiltonian. 

A numerical solution of these equations proceeds as discussed in 
Section III and will not be repeated here. However, several modifica- 
tions of the basic variational problem will be mentioned here. First, 
inequality constraints on the state and/or control variables can be 
introduced and handled by the adjoint technique. This will be the sub- 
ject of a later report. Second, it can happen that some or all of the 
v«'s in (4.29) are all but impossible to eliminate analytically. If 
this happens they can be made part of the convergence process. Assum- 
ing all q of the vg's are to be found, suppose q additional differential 
equations of the form = 0, ..., = 0 are added to the basic set of 

state differential equations. Their solutions are, of course, constants 
whose correct values are those which reduce (4.29) to identities. It is 
easily verified that the Lagrange multipliers and adjoint variables intro- 
duced by these new equations are also constants. Consequently, numerical 
solutions of additional new differential equations are not required. How- 
ever, formally, there are now a total of m + 2q + 1 terminal conditions of 
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which q are of no interest; i.e., q of them are the final values of new 
Lagrange multipliers which do not enter into the problem solution at all. 
Therefore, there are m + q + 1 terminal boundary conditions to be satis- 
fied. This requires q more integrations of the adjoint equations where 
boundary conditions are given by 


P 


(i) _ M 


U) 


a 


Sz 


a 


JZ = 2m - p + 1, . .., 2m - p + 1 + q 
<3=1, ...» 2m + q 


and (3.14a) can again be used to calculate changes in the guessed v^ f s. 
Lastly, the adjoint technique can be used to calculate controls for 
neighboring extremals which can be useful in guidance analysis. 


V. EXAMPLES 


To reduce complexity, the differential equations of motion of the 
flight mechanical problems considered for these examples will be written 
in two dimensions. Further, the reference planet is considered spherical 
and thrust levels are constant. Under these conditions, the equations 
governing the motion of a rocket vehicle operating in a vacuum (i.e., 
exo-atmospheric) can be written 


F GM 

V = £t = — cos a - — 5 COS •a 
x m r 

= f 2 * sin a + v sin -9 (5.1) 


r = f 3 = v cos 9 
m = f 4 = k 


where 


F is the rocket thrust, 

GM is the product of the universal gravitational constant 
and the reference planet mass. 
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v is the velocity, 
m is the instantaneous rocket mass, 
k is a given constant, 

r is the distance of the rocket from the center of the 
reference planet, 

•9 is the angle between the radius vector r* and the velocity 
vector v, 

— ) 

a is the angle between the thrust vector F and the velocity 
vector v*. 

Geometrical relations between these quantities are illustrated below. 



There are a total of 8 variables, namely, r, v, ■9, m, F, GM, k, and a 
of which 3 - F, GM, K - are known constants. Since v, 9, r, and m are 
determined from (5.1), there is one free variable a, the control variable, 
which can be varied to minimize the time of transfer between the given 
initial state and the given terminal state. In the development of the 
Euler-Lagrange equations and the adjoint equations, the following corre- 
spondences for the subscripts and variables are made: 


1 

2 

3 

4 


=> V 

x x => V 

=> 9 

x 2 => 9 

=> r 

x 3 => r 

=> m 

x 4 => m 
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which are read. Subscript 1 is replaced by v, subscript 2 is replaced 
by d, etc. The Hamiltonian (4.24) becomes 


H = ! + q v 


F 

„ GM 



~ F 

/ GM 

l\ 

— cos 
in 

a •- cos 

r 

•a 

+ q 3 

— sm a + 
mv 

l,PP ■ 

■ — v s in 3 

V J 


+ q r (v cos -3) + qjc. 


(5.2) 


The Euler-Lagrange equations (4.26) are 


• ( F . . GM . q. , 1 . 

% “ sin a + 7^ s m *8 + ^ sm -3 J - cos q 


• GM / GM 1 \ 

= " p? sin -3 q v - v cos 3 q^ + v sin -3 q r (5.3) 


q r = - 


2GM 


_ , f ZLri v l 1 \ 

“ C0S 4 % + ■ r) 


2GM 1 \ v sin -3 


r q 3 


Sn 


COS 


a 


m 


<L„ + 


m^v 


sin a q 


•3 


The equation for the control variable (4.27) is 


F F 

- ~ sin a q + — cos aq. 
m n v mv n 3 


= 0 


or 


tan a 


(5.4) 


The terminal boundary conditions (4.28) are 


ft 


(D = 


v - 



0 


ft 




0 


ft 




0 . 


(5.5) 


The terminal boundary conditions (4.29) and (4.30) are 


q 


v 



= 0 


q 





= o 


(5.6) 


q + v 




= 0 



(5.7) 



(5.8) 


21 



Since (5.1) and (5.3) are a total of eight differential equations, ten 
boundary conditions are required. Five of these are given by the assump- 
tion that the initial boundary is known. Apparently, the remaining five 
come from (5.5) - (5.8). Equations (5.5) must be satisfied so that the 
remaining two conditions can be picked from (5.6) - (5.8). Actually, 
equations (5.6) yield no new information. Hence, the five additional 
boundary conditions occur on the terminal boundary and are given by 

(5.5) , (5.7) and (5.8). A solution of the complete differential system 
proceeds by guessing tf and the missing initial boundary values and vary- 
ing them until the integration of (5.1) and (5.3), using (5.4) for the 
control, satisfies (5.5), (5.7) and (5.8) at some tf = tf. 

In the just outlined solution procedure, tf and all four initial 
multiplier values were needed. Two artifices can be used to reduce 
these five initial guesses to three guesses, namely, tf and two of the 
initial multiplier values. This can be done by noting that, first, the 
Euler-Lagrange equations are homogeneous in the multipliers and, second, 
that the first three of the Euler-Lagrange equations are independent of 
the fourth. Homogeneity means that for any solution set {(q v (t), q^(t), 
q r (t), %t(t) {kq v (t), kq^(t), kq r (t), kq m (t) } is also a solution set 

where k / 0 is any real number. Consequently, for the multipliers solv- 
ing the problem, there exists a k such that one of their initial values 
may be fixed at some number N. (N is often given the value 1.) For 
example, A (t Q ) f 0 of a solution set becomes 1 if k = 1/A(t 0 ). Requir- 
ing this condition at the outset reduces the initial guesses by 1. The 
second observation above, along with relaxing the requirement that (5.7) 
be satisfied, means that the initial value of q m (t Q ) may be fixed at 
some convenient value, say 1. The motivation for doing so is that now 
there are three arbitrary parameters, tf, q v (t Q ) and q$(t Q ) , say, to be 
picked (if a solution exists) such that the three boundary conditions 

(5.5) are satisfied. This handling of q m introduces the simplification 
that now equations (3.14a) of the adjoint technique are three in number 
and easy to solve for the required initial corrections. Additionally, 
q m (t) can still be used in computing the integration constant H which 
is an aid in checking the accuracy of the integration algorithm. (Note 
that now H ^ 0.) The multiplier q m could have been eliminated entirely 
from consideration, since by assumption m is a known function of time, 
but this would have introduced an explicit time dependency into the state 
differential equations which is not allowed in the formulation in Sec- 
tion II. (However, t can be introduced explicitly into (2.1) and (3.1) 
with no resultant difficulty in the analysis whatever.) It should be 
remembered that these contrivances for reducing the initial guesses from 
5 to 3 are for convenience and expediency in solving the problem and 
that formally all five should be used to satisfy (5.5), (5.7) and (5.8). 
(It could be said that the three guess solution finds incorrectly scaled 
multipliers while yielding the correct control program and the minimum 
transfer time.) 
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The adjoint equations (4.31) and (4.32) become, after some obvious 
cancellation and rearranging (p* and p9 are ignored for the reasons given 
above) : 


=> x = - 
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(5.9a) 



GM . 0 

X 

/ GM A 


— •g' sm -3 

L r J 

p v - 

" r) V COS 6 


X - . X 

p^ + [v sin 3] p r 


i GM , l\ • a 

q^ + rj cos 3 + q r sm 3 


p q + 

*v 


GM 


q v T5 cos 3 


q^ - “J v sin 3 - q r v cos 3 


P., “ 


2GM . „ 

73 “ sin -8 


* ( 2GM 1 v cos 3 

+ q A ^'7 


(5.9b) 
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\ 


= [cos -3] - [v sin $] p^ . 


(5.9f) 


As mentioned earlier, even for the simple problem here, the adjoint equa- 
tions are fairly lengthy. A solution can proceed by fixing q r (t 0 ), guess- 
ing values of tf, q v (t Q ) and q.s(t 0 ) and integrating (5.1) and (5.3) forward 
in time until t = tf. At this time, (5.5) will normally not be satisfied 
and (5.9) along with (5.1) are integrated backward in time until t = t Q . 

The starting values for the adjoint variables at t = tf are calculated 
from (3.9a) to be 


(l)x - 1 

P v - 

(i)x - n 
" °’ 

( e - o. 

( P l)q - 0, 
r v 9 

- o. 
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(|)X _ 
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- ». 

( # )q = 0 
r 

(5.10) 

(3)X 

p v = o. 

(3) X 

P d = 0 » 

(a)x 

p r - 1. 

( 3 )q - 0 

p — u, 
*v ’ 

^ - 0, 

(3>q n 

p = 0. 



At time t = t Q , the integrated values of the adjoint variables correspond- 
ing to the initial values of (5.10) are substituted into equations (3.14a), 
which become 


Xt +( ^ q "Ot + * dt£ )t * L v ■ v X t • 

c ° c o c f c f 


*0 t +<p ? *.) 
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i +fldt. - -3 - 0. 1 • C 2 

to /t f L Vt f J 


(5.11) 


^ q +<p T X)*. + f dt X f “ [ r - X-j • c ’ 


Where 0 < C x , C 2 , C 3 g 1. 
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The right-hand sides of (5.11) represent a fractional portion of the dif- 
ference between the desired terminal state values and the actual terminal 
state values for this particular integration. Equations (5.11) are solved 
for the increments Sq v (t 0 ), Sq,g(t 0 ), dtf and a forward integration is 
attempted with new initial values 


new 

<lv (t o> 


new 


old 


W 

old 

W 




and a new cut-off time 


new 

t„ 


old 

t f +dt £ . 


(5.12) 


(5.13) 


If the fractional terminal state violations have not been reduced by 
this forward integration, the increments just computed are reduced by 
some factor, say 1/2, (5.12) and (5.13) are recomputed and a new forward 
integration is made. This step is done as many times as required to 
achieve a reduction in the terminal state violations. (If the problem 
has a solution, there will be a reduction after finitely many steps.) 

If equations (5.5) are again not satisfied, a new cycle of backward and 
forward integration is initiated, the whole process being done as many 
times as necessary to achieve convergence. (C^, C 2 , and C 3 are increased 
as the linear predictions become increasingly more accurate.) 

The first numerical example considered is that of transferring a 
rocket vehicle from a 100 nm (185.2 km) altitude circular orbit to a 
300 nm (555.6 km) altitude circular orbit with the earth as the only 
gravitational body. The assumed vehicle and earth model parameters are . 
given in Table 1. 


26 


TABLE 1 


Vehicle 

F = 50 lbf 

Initial Weight = 500 lbf 
Isp = 400 sec 

Earth 

Radius = 6370 km 

GM = 3.98059389 x 10 5 km 3 /s 2 


Practically speaking, the engineering units as typified in Table 1 are 
not the best units in which to solve the problem. In order that the 
various elements of the terms in the adjoint equations be taken account 
of most accurately, it is advantageous to scale the units so that the 
problem variables are of the same order of magnitude. Ordinarily this 
is accomplished by writing the equations in an appropriate nondimens ional 
form. Here, however, the nondimens ional ized form is achieved by properly 
scaling the various input variables. In this way an engineering solution 
can also be obtained by unsealing these input variables. A convenient 
scaling is achieved by adopting the earth 1 s radius as the unit of length 
and requiring GM — 1 in the scaled units. This results in 805.81475 sec 
as the unit of time. The vehicle 1 s initial mass can be used as the unit 
of mass. With these scaling factors, the initial state is 


r = -M i r - = 1.0290738, v = %/l/r - .98577260 

o 6370 o o 


■d = 90°, m = 1. 
o o 


The terminal state is 


i- = = 1.0872214, v, = JT/rZ = .95904947, - 90°, 

f oJ7 0 
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to be maximized; q r (t Q ) and q m (t 0 ) were kept fixed at q r (t Q ) = +1 and 
q m (to) = +1. The initial values of q v and q^ were guessed at q v (t Q ) = +1 
and q$(t Q ) = +1 . The cutoff time was initially guessed as .7 
scaled units. The C f s of equation (5.11) were initially set at .03 and 
were doubled or halved as the iteration proceeded according to whether it 
took more or less than three forward integrations per cycle to reduce the 
terminal boundary violations. The terminal boundary conditions were 
achieved to 8 decimal places within 10 trials. Using an integration time 
step corresponding to 20 sec, the Hamiltonian was constant to 5 decimal 
places using Runge-Kutta 4th order integration formulas. The final cut- 
off time was 1.3718816 scaled units indicating that the initial guess was 
wrong by a factor of .5. Figure 1 summarizes the control history as the 
trials progressed. Figure 2 summarizes the r - v history of the trials. 
The converged values of q v and q^ were q v (t Q ) = .38828166 and q,g(t 0 ) » 
-.54702046. Both figures indicate that the initial guesses were quite 
bad. 


The second numerical example involves this same vehicle in an 
escape mission. The initial boundary is the same and the terminal state 
was chosen to be rf = 2.0175824, v f = 1.6, = 45°, m f maximized. The 

energy of this terminal state is typical of a low energy mission to Mars. 

The initial values of q v and q^ were guessed as +1 each, and the cutoff 
time was guessed as 2.23. The C's were initially set at .015 and sub- 
sequently halved or doubled as before. This time terminal boundary con- 
ditions were achieved to 8 decimal places in 18 trials. The Hamiltonian 
was constant to 5 decimal places using a time step equivalent to 40 seconds. 
The final cutoff time was 4.4778898, while the converged values of q v and 
were q v (t 0 ) = +.23594021 and q^(t Q ) = -.49928473. The control histories 
are summarized in Figure 3 and r - v plots are given in Figure 4. 

A third numerical example, picked from the literature, is an Earth- 
Mars interplanetary transfer mission. Earth and Mars are assumed to be 
in coplanar, circular orbits about the sun. Adopting the astronomical 
unit as the unit of length, and requiring the heliocentric GM to equal 1 
in the scaled system leads to a time unit of 58.13504 days. The vehicle 
parameters are shown in Table II. The initial state is r Q = 1, v Q * 1, 


TABLE II 


Vehicle 

F = .127 lbf 

Initial Weight = 1500 lbf 
Isp = 5700 sec 
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Control Anglo, a (dog) 
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Radius (scaled units) 
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Velocity (km /sec ) 

FI6. 2. RADIUS VERSUS VELOCITY PLOTS FOR TRANSFER 
FROM 100 NAUTICAL MILES (189.2 KM) 

TO 300 NAUTICAL MILES (555.6 KM) CIRCULAR ORBITS 



Control Anglo, 



40 60 80 100 120 140 160 180 200 

Tine (day*) 

CONTROL PROGRAMS FOR HELIOCENTRIC TRANSFER 
FROM EARTH'S ORBIT (ASSUMED CIRCULAR) 

TO MARS' ORBIT (ASSUMED CIRCULAR) 



Radius (scaled units) 
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RADIUS VERSUS VELOCITY PLOTS FOR HELIOCENTRIC TRANSFER 
FROM EARTH'S ORBIT (ASSUMED CIRCULAR) 

TO MARS' ORBIT (ASSUMED CIRCULAR) 



•0 O = 90°, m D = 1. The terminal state is rf = 1.525, v f - .8098, *8f = 90°, 
and mj maximized. Again, q v and were initially set equal to +1 each 
and the cutoff time guessed as 1.66. The C 1 s were initially set at 
.00015 and the integration time step was set equivalent to 5 days. 

Boundary conditions were achieved to 8 decimal places in 26 trials with 
the Hamiltonian being held constant to 4 decimal places. The final cut- 
off time was 3.3194012 with the converged values of q v and q^ being 
q v (t ) * 1.0784030 and q^(t Q ) = -.49498598. Typical control histories 
are shown in Figure 5, and r - v plots are shown in Figure 6. 

The fourth and last numerical example involves launching a two-stage 
vehicle into a low circular earth orbit. The equations of motion are 
written in a two-dimensional earth-fixed coordinate system. The earth is 
assumed spherical and its atmosphere is modeled by an exponential function. 
The drag and lift coefficients of the vehicle are assumed to be constants. 
The differential equations are 


. ^ F GM 0 , 

v = f 1 = — cos CC - — 5 cos *6 + 
x m r 


a) 1 2 r cos fl - ® cos a - - sin a 
m m 


4 = fp s — sin a + 
^ mv 
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-) v sin 4 + 
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— cos a - — sin a 
mv mv 



sin *6 - 



(5.14) 


r = f 3 = v cos -d 
m = f 4 = k. 


where 


D = axial force = 1/2 pv 2 A s 


N = normal force = 1/2 C n pv 2 A s C£ 


p = atmosphere density = p Q e 
u' = (j cos 9 sin A z 


-Q(r-r e ) 
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Cj is the axial force coefficient 
u o 


C n is the normal force coefficient 
u o 

A s is the reference frontal area of the vehicle 

p Q is the atmospheric density at the earth’s surface 

Q is the inverse of the scale height 

10 is the earth’s rotational rate 

r e is the earth’s radius 

cp is the launch latitude 

A z is the launch azimuth 


and the other terms are as defined in (5.1). There are a total of 17 
variables - r, v, ■a, m, F, GM, u', cp, A z , Cd Q , C no , A s , Po, Q, r e , k, a - 
of which 12 (F, GM, w', cp, A z , Cd Q , C no , A s , p Q , Q, r e , k) are known con- 
stants. Equations (5.14) determine r, v, -a, m so that again one free 
variable Ct is available to minimize the transfer time from the given 
initial state to the given terminal state. Using the subscript notation 
defined previously, the Hamiltonian (4.24) becomes 


H = 1 + q ( — cos a - cos •a + to' 2 r cos •a - — cos a - - sin a 
n v \m r^ mm 


+ q. 


F / GM 1 \ N D 
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rnv \r^ ry mv mv 


, 2 

oj r „ „ , 

sin a - 2ur 


+ q r (v cos a) + q Itl (k), 


(5.15) 



The Euler-Lagrange equations (4.26) can be written 
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The control, a, is obtained from (4.27) which becomes 


F F 

— sin a q + — cos a q _ + 

m rav ^-Q 


r A s p 

(Cd o " Cn o ) sin a ' Cn o° ; cos a K“^r j n v 


+ {< c n 0 - c d Q ) cos a - Cn o a sin aj- (jisr) = 0. (5.17) 


Notice that (5.14), (5.16) and (5.17) are identical with (5.1), (5.3) 
and (5.4) except for the bracketed terms which represent the effect of 
the earth’s rotation, the assumed atmosphere and vehicle aerodynamics. 

The terminal boundary conditions are given by (5.5). The adjoint equa- 
tions are much too lengthy to be listed here. However, the partials 
necessary for substitution into their general formulas, (4.10), are 
listed in Appendix II. 

The vehicle, aerodynamic, and launch parameters are given in 
Table III. The formulation so far has not allowed for discontinuities 
in any of the variables, and thus is not directly applicable to a vehicle 
having stages with different characteristics. The vehicle described by 
Table III has discontinuities in thrust and mass at the fixed staging 
times which cause the differential equations involving these two variables 
to also be discontinuous at these times. However, these types of discon- 
tinuities introduce no difficulties since it has been demonstrated in the 
literature that the Lagrange multipliers are continuous across fixed stag- 
ing times when these times are independent of the state. The Hamiltonian 
is discontinuous at these staging times, but the amount is unimportant and 
it still is a constant for each stage. The demonstration of these facts is 
straightforward and starts by defining functions dependent on the state 
variables and time which determine the staging times and the magnitude of 
the discontinuities at these staging times. These functions are adjoined 
to (2.4) with new Lagrange multipliers, and the integral appearing in 
(2.4) is divided into parts over which its arguments are continuous. The 
first variation of the new (2.4) obtained is set equal to zero, and an 
interpretation of the result yields the preceding statements for the 
special case of fixed staging times. The main result of this discussion 
is that the solution technique already discussed holds. Care simply must 
be taken at the staging points that the discontinuous differential equa- 
tions are handled correctly. 
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Vehicle 

1 st stage 

Initial Weight 

1,000,000 lbf 

F 

1,600,000 lbf 

Isp 

300 sec 

Staging Time 

105 sec 

2 nd Stage (Coast) 

F 

0 lbf 

Staging Time 

115 sec 

Weight Drop 

140,000 lbf 

3 rd Stage 

F 

200,000 lbf 

*sp 

420 sec 

Aerodynamic 

A s 25 m 2 

p Q . 13133546 

(kgm/m 3 ) 

Q .13623243 

x 10“ 3 (1/km) 

<*o -5 

<*>o 5 

Launch 

a) .72921157 x 10' 

i 

’ 4 (rad/sec) 

0 28.28° 
A z 72° 





The principal difficulty introduced by adding the atmosphere is the 
solution of (5.17) for the control. It has not been mentioned thus far, 
but the correct value of a computed from (4.27) is the value that maxi- 
mizes the Hamiltonian. This statement can be proven by an application 
of the Pontryagin maximum principle or the Weierstrass E-function test. 
Since the solutions of (5.4) are periodic with period it, at most two 
solutions are possible in the interesting range -n < 0! ^ jt. One of 
these maximizes the Hamiltonian, the other minimizes it. For the vacuum 
flights (5.4) does indeed maximize the Hamiltonian (as long as q r (t Q ) > 0). 
The situation for (5.17) is quite different. First, solutions are not 
periodic with period rt (except in the limit for large a). However, there 
are multiple solutions in the range -it < a § it. The details of the solu- 
tion of (5.17) will be discussed next. 

Equation (5.17) is of the form 


(A + Ba) tan a = C + Da 


(5.18) 


where A, B, C, D are functions of v, m, F, Cd Q > C n Q> A s , qv and *1^* T ^ e 
left-hand side (L.H.S.) of (5.18) is the product of a linear function and 
a trigonometric function. The linear term is called the magnification 
factor since (A + Ba) has the principal effect of increasing tan a in 
absolute value. A graphical representation of the L.H.S. of (5.18) 
divides essentially into four different cases: A 3 0, B > 0 and A < 0, 

B ^ 0, An important subcase occurs when A + Ba = 0 and |a| = n/2 simul- 
taneously. A typical representation for A > 0 and B ^ 0 is shown in 
Figures 7a and 7b. A solution of (5.18) is thus represented by the 
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Figure 7 



intersection of a straight line (the right-hand side (R.H.S.) of (5.18)) 
and the dotted lines (the L.H.S. of (5.18)) in Figure 7 . Depending on 
the values assumed by C and D, there are at most four solutions. Further, 
there are at most two solutions per quadrant. Notice that always there 
are two solutions. The special case A + Ba = 0 and |a| = rt/2 requires 
special treatment. Taking a = -«/2, then A + B(-rt/2) = 0, implies 
A ■ B(ji/2). Therefore, (A + Ba) tan a becomes 


B 



tan a. 


This is an indeterminant of the form 0 • -» as 


Q£ — > *■ 


iff 

2 * 


However, its limit does exist and is evaluated by L* Hospital's rule as 
follows: 



* - B lim sin^ 


a -> - 
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A similar evaluation for 
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n- 
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is continuous everywhere in [-it, 0], Thus, Figure 7a would become Fig- 
ure 8.-. A similar result would hold for Figure 7b. The conclusions 
from Figure 7 still hold. A solution of (5.17) proceeds on the basis of 



Figure 8 



4 


these conclusions. dH(a)/da, i.e., (5.17), is evaluated at increments 
of it/2 in the range [-it, it]. If its sign changes at the end points of 
any of the subintervals, a root exists within this subinterval. If the 
sign of dH/da does not change, the sign of d 2 H (a) /3a 2 is examined. If 
it changes, there are two roots in that sub interval. These two roots 
are placed in appropriate smaller subintervals by evaluating dH/da within 
this subinterval. A modified Newton-Raphson procedure is then used to 
find roots within the intervals in which they have been isolated. The 
modification is as follows. Hie normal Newton-Raphson iteration for a 
root has the form 


a ,= a 

n+1 n 


dH 

W 


The modified form is 


a 


n+l 


dH 

Sol 

“n " KlK2 ’ 

^n 


where Kj, * ±1 and 0 < K 2 § 1. K x and K 2 can be given geometric interpre- 
tations, but it will suffice here to say that K 2 limits the sequence of 
OC 1 s to lie within the interval the root has been isolated within and Kj. 
causes the sign of the derivative to agree with the sign of the finite 
increment slope found in the process of isolating the root. In this 
manner, all the roots are found in the interval -it < a S it. Each is sub- 
stituted into the Hamiltonian and the one maximizing H is picked as the 
control. A typical variation of H(a) is shown in Figure 9., In use, 
the foregoing procedure is very rapid and the three roots indicated in 
Figure 9 would be found to eight significant figures in four or five 
iterations per root. 

The mission for the two-stage vehicle described in Table III is the 
attainment of a circular orbit at 194.6 km altitude. Using the scaling 
given in the first numerical example, the initial state is r Q = 1.0004223, 
v Q =* .020240245, -d Q - 6°, and w 0 = 1. The terminal state is Tf = 1.0305494, 
Vf = .98506657, = 90° and Wf maximized. The initial values of q v and 

q^ were guessed as +.1 and 0, respectively, and the cutoff time was guessed 
as .67. The C’s were initially set at .0001. Terminal boundary conditions 


43 


V 



Figure 9 


were achieved to 8 decimal places in 29 trials with the Hamiltonian being 
held constant to 5 decimal places. A time-step equivalent to 10 seconds 
was used. The final cutoff time was .7906737, and the converged values 
of q v and q^ were q v (t G ) = .074358687 and q^(t Q ) = .0032090726. Control 
histories are given in Figure 10, and r - v plots are shown in Figure 11. 


VI. CONCLUSIONS 


The adjoint method has been shown to be a powerful tool in the solu- 
tion of two-point boundary value problems. In particular, for variational 
problems treated via the Lagrange multiplier technique, solutions are 
freed quite drastically from dependence on the initial multiplier values. 

With respect to the numerical examples of the text, the number of trials 
necessary for convergence in the second and third examples could have 
been easily cut in half by guessing even remotely reasonable cutoff 
times. For the atmospheric example, the number of trials could have been 
reduced considerably by using a smaller time step in the integration 
algorithm since terminal boundary conditions had been satisfied to 2 
decimal places after 15 trials. Generally, the convergence process pro- 
ceeds more rapidly the smaller the integration time step. This is because 
the influence functions are obtained more accurately. The time step required 
for accurate forward integration is generally not the same as that for accur- 
ate backward integration. There is a trade-off between how accurately the 
influence functions need to be obtained to solve the problem and how many 
trials are required. Usually, the time-step for a backward integration 
needs to be smaller than for a forward integration for equivalent degrees 
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Control Angle, a (deg) 
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FIG. 10. CONTROL PROGRAMS FOR ASCENT 

THROUGH AN EXPONENTIAL ATMOSPHERE 

105 NAUTICAL MILE (194.6 KM) CIRCULAR ORBIT 




Radius (scaled units) 
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FIG. 11. RADIUS VERSUS VELOCITY PLOTS FOR ASCENT 

THROUGH AN EXPONENTIAL ATMOSPHERE 
TO A 105 NAUTICAL MILE (194.6 KM) CIRCULAR ORBIT 




of accuracy. Two disadvantages of the adjoint technique are the slight 
programming problem caused by the backward integration and the time 
penalty incurred, because the state equations must be integrated back- 
ward as well as forward (or their forward values saved on tape and then 
interpolated during the backward integration). An almost one- third com- 
putational time saving can be achieved in launch and ascent problems by 
using velocity as the cut-off criteria. Velocity is an acceptable function 
for this purpose since it is usually a monotonically increasing function 
of time in this type of problem. 
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APPENDIX I 


In the ordinary differential calculus, the notation dy denotes an 
infinitesimal, i.e., a change in the variable which can be made as small 
as desired. In the calculus of variations, Sy denotes an infinitesimal 
change in the function y; i.e., Sy is a variation of y(x). If the original 
function is denoted y(x) and the changed function denoted y(x), then 
y(x) * y(x) + Sy . The change (variation) can be made more explicit by 
saying y(x) = y(x) + k0(x), where k is an arbitrary constant and 0(x) is 
independent of k and continuous over the same range that y(x) is continuous. 
For example, the straight line y(x) = x, 0 ^ x ^ 1, can be deformed into 
the parabola y(x) = x 2 , 0 ^ x £ 1, by defining 0(x) * x 2 - x. Then, 
y(x) * x + k(x 2 - x). Thus, when k = 0 the original function is retained. 
When k = 1, the varied function is the parabola. This little example 
illustrates how a given function can be deformed into another given 
function. Ideas are fixed more firmly with the following definitions: 

Def. 1 Let y(x) and y(x) be uniformly continuous functions in the 
interval xq ^ x ^ x x . Then y(x) is said to lie in a strong 
neighborhood N € of y(x) if and only if for every e > 0, 

|y(x) - y(x) | S e. 

Def. 2 Let y(x) and y(x) be differentiable, uniformly continuous 
functions in the interval x Q ^ x ^ x x . Then y(x) is said 
to lie in a weak neighborhood N e of y(x) if and only if for 
every e > 0, |y(x) - y(x) I ^ € and 

1^ yOO - £ y(x) I * €. 


Strong variations are associated with Def. 1 and weak variations are 
associated with Def. 2. Futhermore, variations are (1) special if the 
independent variable is not varied and (2) general if the independent 
variable is varied. Almost exclusively special, weak variations (although 
very often this is not stated explicitly) are used in deriving the neces- 
sary conditions for extremals in the variational calculus. 


Now consider a function of more than one argument, say £(x, y x , ., 


y m > 


y{ 3 


y m ), where yf = — yi- 


J± . Consider the effect of replacing 

each of the arguments by its varied value; i.e., y^ is replaced by 
yi + MiCx), yj is replaced by y[ + k0[(x), etc. If £(x, y t , .... y m , 
y{ $ •••> y m ) has a series representation, then 


£(x, y£ + k0 i} y[ + k0[) = £(x, y*, y[) + 6^ + j 5 2 £ + ... 
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where &£, 5 2 £, . .., S n £ are the 1 st , 2 nc *, and n*"* 1 variations of £ and 

symbolically 



Substituting successively y-L , . .., y m , yj_ , y m for £ yields 

5yi = b 2 ?! =0 ... fjHy! = 0 

5 y m = k $n &% m = 0 ... S n y m = 0 

8y 1 = =0 ... 8 n y t = 0 

• » ♦ 

& y m * k #n “ 0 • • • = 0 


from which follows easily that Sy* = k£h = — (kj^) = ~ (5y.). 

J J at J at J 


When the end point of a variational problem is not fixed but instead 
is allowed a general variation, the total change in the end point value 
is found as follows. Referring to the drawing, y(t) is the varied curve 
and y(t) is the original curve; y(t) must be extended over the interval 
At. Then, since y(t) and y(t) are assumed differentiable. 


50 


y 

t 



t 


y(t Q + At) = y(t D ) + y(t Q ) At + e x At (1-1) 


where lim = 0. 

At -» 0 


Since Def. 2 holds, 


|y(t 0 ) ■ y(*o> I < e 


or 


y(t 0 ) * y(t 0 ) + € 2 


( 1 - 2 ) 


where lim e 2 = 0. 

Substituting (1-2) into (1-1) yields 

y(t 0 + At) = y(t Q ) + y(t Q ) At + (e^. + e 2 ) At. 

Adding and subtracting y(t Q ) yields 

y(t Q + At) = y(t Q ) - y(t Q ) + y(t D ) + f(t Q ) At + + e s ) At 

or 

y(t Q + At) - y(t D ) = 8y(t 0 ) + y(t Q ) At + (e x + e 2 ) At, (1-3) 

where Sy(t c ) = y(t Q ) - y(t D ) by definition. (1-3) shows that the principal 
part of the total change in the end point is 

Ay = y(t D + At) - y(t D ) = 5y(t 0 ) + y(t D ) At. 
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APPENDIX II 


The differential equations for the influence functions are obtained 
by substituting the following partials into (4.22) and (4.23): 


= (Cd 0 cos a + Cnoa sin a) 



bv <31 


sin -6 - a)* r sin •6. 


cos -9 + (Cd Q cos CC + CiiqQ! sin Q!) 2^““J ^ 


i* 2 cos d. 


dv F cos 0 ! 

"5m 3 


3T 4 + (Cd Q COS a + Cn o a sin a) 


dv F . 

S " “ 3 sin 


« + (< c d 0 - Cn o ) sina- C^a cos a) (^T“) 



^ * - ^“2 sin a + sin *8 + ^ sin *3^ -J^(Cd 0 sin a - C^a cos a) 

w' 2 r s in -fll 

• V } 


/ GM l\ - .. 0 ) ,2 r cos •a 

5» “ IT 2 ^ ’ 7 ; < v 008 *> 7 

^ - 7) (^) ' + < V cos a - c o .in a) 
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d*3 F . 


= - sin a - (-Cd Q sin a + C no a cos a) (^ 2tn 2 J m 


^ = ““ cos a + ^(Cn 0 " c d D ) cos a “ c n 0 a sin a ^) 2m ) 


= COS ■d 


= - v sin -fl 


dr _ 0 

S" 0 


U - - 5 cos a a - — sin a q* + ( (C do - 2Cn o ) cos a 


+ Cn o a sin OCj 2my J % + [v^o “ 2 Cn o ) Sin a 


- v cos a ) Cit) v 
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c> g H F 

= " PP COS « - 


(c^ - c do ) sin a + c nr oc cos a 


•)C J 


a 


% 


(Cn 0 - Cd o ) cos a - C n a sin a 


^ A P A u' 2 r ! 

c) (- 5 — ) + - - ■ ' 

J \2m J V 


sin $ 
V. 


B 2 H _ /’2F . „ . 2GM . 0 , 

- w* sin a + pp sin V % ■ 


(Cd 0 cos a + Cn 0 Q: sin a) 


<£)] 


2co T r sin -6 


P v 


d 2 H 

c^vd-9 


= _ + I 

\r^v^ 


cos « q - sin « q + £° S * q 

•3 r 


c» 2 H _ f 2GM , l\ /sin 3 




r J 


(C dQ sin CC - Cn Q a cos 


(Cd- COS a + C n a sin a) 

w O 


(«£)] 


% 


A . to 1 2 sin -a 

°° v~y + —p — _ 


d 2 H _ n 

3355“ 0 


S 2 h 


(Cn Q - Cd 0 ) sin a + C^a cos a 


,n O L 





2m 


«v 


( c d Q - Cn o ) cos 


\ /QvA Pn- 

a + c no a sin ctj y ■ — ■ • ) 


‘3* 
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d 2 H GM , f GM l\ . p 

FF = P cos « --Jvsiniq^.vcosflq^io' r cos * 


, co* r sin ^ 

+ v V 


d 2 H = _ 2GM 

d’Sdr " r 


sin ■O 


2GM 1 \ v cos -9 


% ■ iT 2 ^ ■ 7/ 7 - w sin ^ q v - 


CO* cos $ 


l r 


Y 2 H _ 6GM a . / 6GM 2 , 

FP - " pr- cos 6 - -p) v sin 3 - 


(Cd Q cos a 


+ Cn o a: sin a) 


z v z k 


2m 


(Cd„ sin a - Cn a cos a) , «_ 

o o \ 2m 


2 vA_P 


•Yl 
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